package edu.uaskl.jqf.tools.matrix.transformations;

import edu.uaskl.jqf.model.Complex;
import edu.uaskl.jqf.model.matrix.MatrixComplex;
import edu.uaskl.jqf.model.matrix.MatrixTransformationHessenberg;
import edu.uaskl.jqf.tools.ArrayTools;
import edu.uaskl.jqf.tools.MathTools;

public enum Hessenberg {
    /** Singleton */
    INSTANCE;

    /**
     * In linear algebra, a Hessenberg matrix is a special kind of square matrix, one that is "almost" triangular.<br>
     * To be exact, an upper Hessenberg matrix has zero entries below the first subdiagonal, and a lower Hessenberg matrix has zero
     * entries above the first superdiagonal.<br>
     * <br>
     * More information and source: <a href="http://en.wikipedia.org/w/index.php?title=Hessenberg_matrix&oldid=540429194">http://en.wikipedia.org/wiki/Hessenberg_matrix</a>
     * 
     */
    public MatrixTransformationHessenberg transform(final MatrixComplex matrixComplex) {
        if (!matrixComplex.isQuadratic())
            throw new IllegalArgumentException("Matrix is not quadratic: " + matrixComplex.getDimensionAsString());

        if (matrixComplex.getRows() == 2) // then it must be 2x2, because it is quadratic
            return new MatrixTransformationHessenberg(matrixComplex, MathTools.intArrayFilledwithIndicesOfSize(2));

        // based on http://www.csee.umbc.edu/~squire/cs455_l13.html by Jon Squire
        final Complex A[][] = matrixComplex.getValuesCopy();
        final int size = A.length; // checked before call
        final int transformationVector[] = MathTools.intArrayFilledwithIndicesOfSize(size);

        // main reduction loop
        for (int currentColumn = 1; currentColumn < (size - 1); currentColumn++) { // over all columns, excluding the first and last
            int biggestSumAbsIndex = currentColumn;
            Complex biggestSumAbs = Complex.get00();
            for (int j = currentColumn; j < size; j++) { // for the current column, get the field with the biggest sumabs (start with diag.)
                final Complex currentField = A[j][currentColumn - 1];
                if (currentField.getSumAbs() > biggestSumAbs.getSumAbs()) {
                    biggestSumAbs = currentField;
                    biggestSumAbsIndex = j;
                }
            }

            if (biggestSumAbsIndex != currentColumn) { // we have a candidate, rowcol row column interchange of H
                ArrayTools.swapArrayIndex1Index2(transformationVector, currentColumn, biggestSumAbsIndex); // swap transformation index
                ArrayTools.swapArrayIndex1Index2(A, biggestSumAbsIndex, currentColumn); // swap rows

                // swap values from all columns
                for (int currentRow = 0; currentRow < size; currentRow++)
                    ArrayTools.swapArrayIndex1Index2(A[currentRow], biggestSumAbsIndex, currentColumn);

                // change fields
                for (int currentRow = currentColumn + 1; currentRow < size; currentRow++) {
                    Complex y = A[currentRow][currentColumn - 1];
                    if (y.getSumAbs() > 0.0) {
                        y = y.divide(biggestSumAbs);
                        A[currentRow][currentColumn - 1] = y;
                        for (int j = currentColumn; j < size; j++)
                            A[currentRow][j] = A[currentRow][j].subtract(y.multiply(A[currentColumn][j]));
                        for (int j = 0; j < size; j++)
                            A[j][currentColumn] = A[j][currentColumn].add(y.multiply(A[j][currentRow]));
                    }
                    A[currentRow][currentColumn - 1] = Complex.get00(); // just cleanup
                }
            }
        } // end main reduction loop

        return new MatrixTransformationHessenberg(MatrixComplex.createMatrixBasedOnArray2d(A), transformationVector);
    }
}
